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Abstract 

A 3-D space-time CE/SE Navier-Stokes solver 
using an unstructured hexahedral grid is de- 
scribed and applied to a circular jet screech 
noise computation. The present numerical re- 
sults for an underexpanded jet, corresponding 
to a fully expanded Mach number of 1.42, cap- 
ture the dominant and nonaxisymmetric ‘B’ 
screech mode and are generally in good agree- 
ment with existing experiments. 

1 Introduction 

Jet noise is one of the most important topics in com- 
putational aeroacoustics. Many of its aspects are of 
primary practical importance and the associated com- 
plicated physical phenomena are the topic of many ex- 
perimental and theoretical investigations. Refs. [ 1 — 4] 
provide a comprehensive discussion and further refer- 
ences. An over/under-expanded supersonic jet emits 
mixing noise, broadband shock-associated noise, as well 
as screech tones under certain conditions. The mixing 
noise is directly associated with large-scale structures, 
or instability waves, in the jet shear layer; whereas, the 
broadband shock-associated noise and screech tones are 
associated with the interaction of these waves with the 
shock-cell structure in the jet core. The screech tones 
arise due to a feedback loop, i.e., some of the acous- 
tic waves generated by the wave/shock-cell interaction 
propagate upstream and regenerate the instability waves 
at, or in the vicinity of, the nozzle lip. The feedback loop 
leading to screech tones is sensitive to small changes in 
the system conditions, and the understanding of the phe- 
nomena is to date based mostly on experimental obser- 
vations [5-8]. 
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Although substantial progress in numerical computa- 
tion of jet mixing noise has been made, reliable direct 
numerical simulation of jet screech noise has up to quite 
recently not been feasible. Shen and Tam [9, 10] ob- 
tained excellent results in direct numerical simulations 
of screech for circular jets using the well-known DRP 
scheme. In their 3-D computation [10], a spectral method 
was adopted in the azimuthal direction, and by using only 
a limited number of spectral functions substantial sav- 
ings of computer memory and CPU time was achieved, 
without deterioration of accuracy. Other recent compu- 
tational work includes [11], 

The present authors [12-14] successfully computed 
axisymmetric near-held screech noise for round jets us- 
ing the recent CE/SE (space-time conservation-element 
and solution-element) method utilizing (unstructured) 
triangulated grids. Because of the implementation (based 
on flux balance) of the non-reflecting boundary con- 
ditions (NRBC), a much smaller near field computa- 
tional domain can be used with this method. However, 
when the CE/SE method is applied to a 3-D rectangu- 
lar jet screech noise computation [ 15], a major challenge 
emerges. Namely, the number of unstructured tetrahedral 
cells required to achieve a certain resolution could be as 
high as tens of millions, which is currently somewhat be- 
yond the capability of common parallel computers such 
as Linux PC clusters. An alternative is to employ an un- 
structured hexahedral grid. The number of cells may then 
be greatly reduced without much loss of resolution. The 
3-D CE/SE Navier-Stokes (N-S) solver needs to be mod- 
ified to accommodate the unstructured hexahedral grid, 
however. 

A general description of the CE/SE Euler method with 
a hexahedral grid was first given by Zhang et al [16]. In 
the present paper, the current 3-D CE/SE N-S solver us- 
ing an unstructured hexahedral grid is described in §2. 
The parallel computation implementation of the scheme 
is outlined in §3. A 3-D circular jet screech problem 
is described in §4. The numerical results are presented 
and compared with the available experimental data in §5. 
Concluding remarks are presented in §6. 
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2 The Numerical Scheme 

As stated in §1, for fully 3-D computations of com- 
plicated problems such as the jet noise problem, a large 
number of computational cells are needed to provide suf- 
ficient numerical resolution. This imposes severe mem- 
ory and speed requirements on the computer system. Of 
course, the choice of unstructured grid type to be used in 
the computation greatly affects these demands. It turns 
out that the use of an unstructured hexahedral grid in the 
CE/SE scheme rather than the heretofore standard tetra- 
hedral one leads to significant reductions of the memory 
and CPU time requirements. Its drawback, of course, is 
a modest increase in dissipation in the scheme. Zhang et 
al [16] first presented a 3-D CE/SE Euler method using 
hexahedral cells instead of tetrahedral ones and provided 
a detailed description of the implementation. 

In general, the CE/SE method [17, 18] systematically 
solves a set of integral equations derived directly from 
the physical conservation laws and naturally captures 
shocks and other discontinuities in the flow. In order to 
have a compact cell stencil, both conservative variables 
and their derivatives are computed simultaneously as un- 
knowns. A brief sketch of the 3-D CE/SE scheme with 
hexahedral cells is given below. 

2.1 Conservation Form of the 3-D Unsteady 
Compressible Navier-Stokes Equations 

Consider a dimensionless conservation form of the un- 
steady 3-D Navier-Stokes equations of a perfect gas. Let 
p, u, v, w, p, and 7 be the density, streamwise trans- 
verse and spanwise velocities, static pressure, and con- 
stant specific heat ratio, respectively. The 3-D Navier- 
Stokes equations then can be written in the following 
vector form: 

Ut + F x + Gy + H z = 0 , ( 1 ) 

where x, y, and z are the streamwise, transverse, and 
spanwise coordinates, and t is time. The five components 
of the conservative flow variable vector U are given by 

Ux = p, U 2 = pu , U 3 = pv, U4 = pw , 

U 5 =p/( 7 — 1) + p ('« 2 + v 2 + w 2 )j 2 . 

The flux vectors in the x, y, and z directions, F, G, and 
H, respectively, are further split into inviscid and vis- 
cous fluxes, 

F = Fi- F v , G = Gi G v , H = Hi — H v , 

where the inviscid fluxes are 

F n = U 2 , 

Fi 2 = ( 7 — 1 )^ 5 + [(3 - !)U 2 - (7 - 1 )(U 2 + U 2 )} /2 Ux, 
Ft 3 = U 2 U 3 /Ux, Fi 4 = U 2 U 4 /Ux, 


Fib = lU 2 Ub/Ux - (7 - l)U 2 [U 2 + Ut + Ui\ /2 7/ 2 ; 

Ga = U 3 , G i2 = U 2 U 3 /Ux, 

G l3 = ( 7 -l)[/ 5 +[(3 - 7 )Ut - (7 - 1 )(U* + U 2 )} /2Ux, 
G i4 = U 3 U 4 /Ux, 

Ga = lU 3 U 5 /Ux - (7-1 )U 3 [U 2 + U 2 + C/ 4 2 )] /2 U 2 -, 
Ha = U 4 , H l2 = U 2 U 4 /Ux, H l3 = U 3 U 4 /Ux, 

Hu = ( 7 — l)t/5+ [(3 - 7 )Ut - (7 - 1 )(U 2 + U 2 )} /2 Ux, 
Ha = lUxUb/Ux - ( 7 - 1 ) 7/4 [Ut + U 2 + Ut)\ /2U 2 ; 

and the viscous fluxes are 

F v x = 0, F v 2 = p,(2 u x - |A), 

F v 3 = n{uy + v x ), F v 4 = p(u z + w x ), 

F v 5 = p[2 uu x + v{u y + v x ) + w(w x + u z ) - |uA+ 

7 d 7/4 u 2 +v 2 

2 

G v x = 0 , G v 2 = p(v x + Uy), 

G v 3 = p(2 Vy - I A), 

G v 4 = p{v z + Wy), 

G v 5 = p[2vv y + u(u y + v x ) + w(w y + v x )w — |f A+ 

7 d 7/4 u 2 + u 2 

P^^lh 2 

H v x — 0, H v 2 piprix + tiz)i 

Hy 3 + tUy), Hy 4 p(2‘W z 3 A), 

Hyb = p[2ww z + u{u z + w x ) + v(w y + v x ) - §uA+ 

7 d 7/4 u 2 + u 2 

PrdyWx 2 

where u, v, w, u x , u y , u z , etc. are the flow velocities 
and their spatial derivatives, which can be expressed in 
terms of the conservative variables Ux, U 2 , etc. along 
with their gradients. Pr = 0.72 is the Prandtl number, /; 
the kinematic viscosity, and 

A = U X + Vy + w z 

is the velocity divergence. 

By considering (x, y , z, t) as coordinates of a four- 
dimensional Euclidean space, E 4 , and using Gauss’ di- 
vergence theorem, it follows that Eq. (1) is equivalent to 
the following integral conservation law: 

<f I m dS = 0, m = 1,2, 3, 4,5, (2) 

Js(V ) 

where S(V) denotes the surface around a volume V in 
E 4 and I m — ( F r n , G rn , H m , U rn ] . 
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Figure 1 : A hexahedral cell with its six neighboring av- 
erage cell centers (ACs). For clarity, only three of them, 

A, B , C are shown. Together, they form a polyhedron of 
24 faces or the base of the CV. 

2.2 Unstructured Hexahedral Grid 

Fig. 1 illustrates a hexahedron cell 1234 — 5678 in the 
unstructured grid. Each of the six quadrilateral surfaces 
(e.g. 1265, 2376, etc.) is associated with a neighbor- 
ing hexahedral cell. Note that the four nodes of such a 
quadrilateral may not lie on the same plane. Assume A, 

B , C,..., are the average centers (ACs) of the respective 
neighboring cells. The coordinates of the average center 
of a hexahedron are the simple average of the coordinates 
of its eight vertices. The average center may not neces- 
sarily coincide with the geometrical (gravity) center. The 
eight vertices of each hexahedron cell and the six average 
centers ( A , B , C..., only three are shown in Fig. 1) of the 
neighboring cells form a polyhedron of 24 faces, which, 
when incorporating the time t direction, is the finite vol- 
ume V in E,\ in Eq. (2) or the control volume CV. 

2.3 Compact Updating 

For any explicit time-marching scheme, flow data at the 
neighboring nodes at the previous time step are required 
to update the flow data at the current node to the present 
time level. In a finite difference method, a difference 
scheme is utilized to do the updating. In most finite vol- 
ume schemes, the computational cell (e.g., 1234-5678 in 
Fig. 1) is the control volume (CV). First, the fluxes on 
the (hyper) surfaces S(V) (e.g., 1265 in Fig. 1) need 
to be updated. Conventionally, flow data is extrapo- 
lated from several neighboring nodes to the center of 
the hyper-surface. In the CE/SE algorithm, a compact 
updating is employed [17, 18]. The compact updating 
achieves high resolution by using a cell stencil consist- 
ing only the immediate neighboring cells. For exam- 
ple, the CV is expanded from the hyper-hexahedral cell 
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1234 — 5678 to a hyper-polyhedron of 24 faces (Fig. 1), 
with all the six neighboring cell centers (A, B , C...) in- 
cluded as vertices. When applying Eq. (2) to this hyper- 
polyhedron in E±, the flux associated with the face 1265 
is now replaced by the fluxes associated with the four 
triangles AA12, AA26, AA65 and AA51 for the new 
CV. The advantage of this procedure is that flow data 
at the cell center A is known and since A lies in each of 
these triangles, no extrapolation through the interior of 
the CV is needed when evaluating the fluxes associated 
with these triangles. The updating is thus conservative. 
In the CE/SE scheme, in order to achieve higher order 
accuracy, flow data are extrapolated by Taylor expan- 
sion from A to the triangle centers along the triangle sur- 
face plane. Consequently, not only the conservative flow 
variables U but their spatial derivatives also are consid- 
ered as the unknowns (totally 20 scalar unknowns) Note 
that Ut is obtained by evaluating Eq. (1). The hyper- 
24-face-polyhedron is the control volume where Eq. (2) 
is applied for conservative updating. Any of these 24 
surfaces (segments of solution elements or SEs) is as- 
sociated with one of the six neighboring average cen- 
ters (ACs, e.g. A, B,C, ...), where the solution U and 
its gradients are already given at the previous time step 
n. Then both the inviscid and viscous fluxes on these 
24 faces (hyperplane surfaces in E,\ space) can be calcu- 
lated by evaluating the flux functions F, G and H at the 
face geometrical centers through Taylor expansion. Con- 
sequently, the unknown U n+1 at the new time level n+1 
is evaluated at the geometrical center of the polyhedron 
from the divergence theorem (2) in E± space. Note that 
the ‘geometrical center’ or centroid of the polyhedron is 
in general different from its average center (AC). After 
evaluating the new gradients (see next subsection) U x , 
U y , U z , it takes only one more step of Taylor expansion 
to extrapolate U n+1 to the average center (AC). 

The present compact updating avoids the uncertainty 
of dimensional-splitting and extrapolation, and hence 
yields better accuracy. In addition, no Riemann solver is 
needed at these hyper-surfaces. Incorporation of the gra- 
dients of U further enhances the accuracy of the scheme 
and makes it possible to use a compact cell stencil. 

After all the U n+1 at the polyhedron centroids are up- 
dated, the next step is to calculate their corresponding 
spatial gradients. 

2.4 Evaluation of Spatial Gradients 

As in the standard CE/SE procedure [17], unknown spa- 
tial derivatives/gradients U x , U y , and U z are evaluated 
using the weighted average technique or an extended van 
Albada limiter [19]. The six neighboring polyhedral cen- 
ters A, B , C. I). E, F around a hexahedron cell (Fig. 2) 
form an octahedron. The vertices of each of the trian- 
gular face of the octahedron and the polyhedral center O 



Figure 2: An octahedron formed by the six neighboring 
polyhedron centroids A , B 1 C, D, E, F, with the current 
polyhedron cell centroid O in the middle. With O always 
as a vertex, eight tetrahedra are formed and hence eight 
sets of gradients exist. 


associated with the current hexahedral cell form a tetra- 
hedron (e.g. ABEO in Fig. 2). The U n+1 are already 
calculated at all the vertices. For each tetrahedron, a set 
ofC/" +1 , Uy +1 , and C /? +1 can be directly calculated by 
solving a linear equation system. Totally, there are eight 
different tetrahedra, and so also eight sets of the spatial 
derivative data. The final values of the spatial derivatives 
are obtained by applying the extended van Albada limiter 
(weighted average) to these eight sets of data. Let U x ^\ 
U y ®, U z ^\ i = 1 , 2 ,..., 8 be respectively the gradients 
from the /th set data. Let 

n = \{u x ^f + ( u y (i) f + (t/, (i) ) 2 r, 

where a is any real index number. Let t, — t\ ~ 1 ; if 
Ti = 0 , ti = 0 and that 


t — ti + (2 + £3 + + t$, 


then, with the van Albada limiter. 


TT S uj% TT £ Uy^ti TJ EE7,«ii 
U x = : ,Uy = y - ,U Z = 


t 


t 


t 


As is well-known, the van Albada limiter is less dif- 
fusive. In practice, it was found that if the power in- 
dex number a is chosen to be slightly negative, say, 
a = — 0 . 2 , the numerical dissipation resulting from av- 
eraging could be further reduced. 

Once the new gradients are evaluated, U n+1 at the AC 
of the hexahedron is obtained by Taylor expansion from 
the centroid of the polyhedron. One marching step is thus 
completed. 



Figure 3: 2-D quadrilateral mesh in a cross section in a 
3-D hexahedral grid 

3 Hexahedral Grid Generation and Parallel 
Computation 

In order to substantially reduce the required number 
of cells and to increase numerical stability in the compu- 
tation, a hexahedral grid is adopted for the 3-D CE/SE 
Navier-Stokes (N-S) solver. For the simple geometry of 
a circular jet, with the x axis assumed in the stream di- 
rection, the hexahedral grid can be generated as follows: 

1 . a 2-D quadrilateral unstructured mesh is generated 
in a circular domain on the y — z plane (Fig. 3); 

2. the 2-D mesh is translated step by step in the x di- 
rection to form the 3-D hexahedral mesh. At loca- 
tions that are occupied by the jet nozzle body, no 
cells are generated. 

In the present 3-D screech noise computation, employ- 
ing a hexahedral grid instead of a tetrahedral one helps to 
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Figure 4: Schematic diagram of parallel computation. 



jet 
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Figure 5: A typical partition of the computational do- 
main by METIS, with different shading or color indicat- 
ing the subdomains. 


reduce the number of cells from tens of millions to a few 
millions. Still, the number of computational cells is very 
large. Due to the large number of computational cells 
(about 3.67 millions in this case), parallel computation 
with multiple processors becomes necessary in view of 
computation turn-around time and memory size. 

The parallelization procedure is similar to the one de- 
scribed in [15] and is sketched in Fig. 4: 

1 . the unstructured hexahedral grid is generated for the 
entire computational domain; 

2. the domain is decomposed into subdomains accord- 
ing to the assigned number of processes (usually 
one-to-one with CPU’s, here 60 are used), using the 
METIS code. METIS is an efficient mesh partition- 
ing code and is freely available from the University 
of Minnesota [20] — for example. Fig. 5 illustrates 
a typical partition of the computational domain for 
the current circular jet noise problem; 

3. the N-S flow solver is modified to use MPI library 
calls and applied to each subdomain to conduct 
computations, with neighboring domains exchang- 
ing pertinent results. 

MPI, or message passing interface, is an interprocessor 
communication protocol standard. The software library 
is prepared by the Argonne National Laboratory [21]. 

4 The 3-D Jet Screech Noise Problem 

The circular jet in 3-D space is sketched in Fig. 6. The 
flow at the nozzle exit is choked, i.e., the nozzle exit 


Mach number, M e , is unity, and the ambient air is sta- 
tionary. The case of jet Mach number Mj = 1.42 is 
considered. These conditions correspond to the experi- 
mental conditions of Panda [5, 6]. In these experiments, 
it was shown that for Mj = 1.42, the jet noise field ex- 
hibits truly 3-D phenomena, e.g., the flapping ‘B’ mode. 
Hence, a 3-D Navier-Stokes solver described above is re- 
quired. 

In the investigation, our attention is focussed on the 
near field of the nozzle since this is the source region 
of the noise. The inner diameter, D, of the jet nozzle 
is chosen as the length scale. The density, po , speed of 
sound, ao, and temperature To in the ambient flow are 
taken as scales for the dependent flow variables. 

In order to clearly display the upstream propagating 
screech waves, the computational domain was extended 
2D upstream of the nozzle exit. The full computational 
domain is a circular cylinder of 10 D axial length and 
7.5 D radius. At the nozzle exit, the inflow plane is re- 
cessed by two cells so as not to numerically restrict or 
influence the feed-back loop. A straight nozzle lip of 
0.25 D in thickness is adopted in the computations. Note 
that in the experimental setup of Panda [5, 6], the actual 
nozzle diameter D = 25.4 mm and the nozzle exit lip is 
beveled. 

The unstructured hexahedral grid currently used is 
generated as described in §2. The hexahedral cells are 
non-uniform since good grid resolution is needed in the 
jet shear layer and in the screech feedback loop paths. 
Cell numbers in the x and radial directions are typically 
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Figure 6: Geometry of the computational domain: the 
inner diameter of the nozzle is D ; the nozzle extends 
into the computational domain by S = 2D, the domain 
length L = 8 D and the radius is R 5.5 D (sponge 

zones are excluded) 

400 and 160 respectively. In the azimuthal direction, 
there are 64 cells, corresponding to totally about 3.67 
million hexahedral cells. The last 6 cells in the stream- 
wise and radial directions have exponentially growing 
size, respectively, and serve as sponge zones to substan- 
tially eliminate any small remaining numerical reflec- 
tion from the outflow and circumferential non-reflecting 
boundaries. 

4.1 Initial Conditions 

Initially, the flow of the entire domain is set at the am- 
bient flow conditions (using nondimensional variables), 
i.e., 

Po = 1, Po = 

7 

uo = 0 , v 0 = 0 , w 0 = 0 . 

4.2 Boundary Conditions 

At the inlet boundary, the conservative flow variables and 
their spatial derivatives are specified to be the same as the 
ambient flow, except at the nozzle exit, where an elevated 
pressure is imposed, i.e., the jet is under-expanded, as in 
the physical experiments. By using the ideal gas isen- 
tropic relations, it follows that the nondimensional flow 
variables at the nozzle exit, with M e = 1, are given by 

_ 7(7 + 1 )Pe 
Pe 2 T r 

2+( 7 -1)M j 2 1^ t 
7 + 1 

( 2 T r \ 

U e = , V e = 0, W e = 0, 

V7 + 1/ 


where T r is the reservoir (plenum) temperature. We will 
also follow the experimental cold-flow condition where 
the reservoir temperature equals the ambient one, i.e., 

T r = 1. 

At the circumferential and outflow boundaries, the 
Type II and Type I CE/SE non-reflecting boundary con- 
ditions as described in the next subsection are imposed, 
respectively. The no-slip boundary condition is applied 
on all the nozzle walls. 

4.3 Non-Reflecting Boundary Conditions 

In the CE/SE scheme, non-reflecting boundary condi- 
tions (NRBC) can be easily constructed based on plane- 
wave propagation theory for hyperbolic conservation 
laws [22], There are various implementations of the 
non-reflecting boundary condition (NRBC) and in gen- 
eral they have proven to be well suited for aeroacoustic 
problems [14, 23]. The following 3-D NRBCs are em- 
ployed in this paper. 

For a ghost grid node [j, n) lying at the outer radius of 
the domain the non-reflective boundary condition (Type 
II) requires that 

(U X )] = ( Uy )7 = (U Z )] = 0, 

while t/" is kept fixed at the initial steady boundary 
value. At the downstream boundary, where there are 
substantial gradients in the radial direction, the non- 
reflective boundary condition (Type I) requires that 

(u x ); = 0, 

while [/”, ( U y )" and (t 7 Z )" are now defined by simple 
extrapolation from the nearest interior node j', i.e., 

U7 = u]„ 

(u y )? = {u y )$, (U Z )? = (U Z )],. 

As will be observed later, these NRBCs, when combined 
with the above buffer/sponge zones, are robust enough 
to allow a near field computation without disturbing or 
distorting the flow and acoustic fields. 

5 Numerical Results 

In this section, 3-D numerical results for the under- 
expanded circular jet described above are presented and 
compared to experimental results [5, 6]. Computations 
are conducted for the jet Mach number Mj = 1.42. At 
this moderate jet Mach number, the dominant unsteady 
motion in the experiments, see [5], is truly three dimen- 
sional. With an appropriate time step size, a large num- 
ber (510,000) of time steps are performed in order to 
achieve sufficient accuracy in the Fourier analysis of time 
series data. Note that no harmonic forcing is imposed 
in the numerical simulation. The initial impact of the 
boundary condition at the nozzle exit stimulates the jet 
shear layer and triggers the feedback loop that generates 
the (then) self-sustained screech waves. 


1 

Pe = ~ 

7 
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Figure 7: Time-averaged experimental (top) and instan- 
taneous numerical (bottom) Schlieren pictures, showing 
shock cell structure; Mj = 1.42. 


Figure 8: Pressure iso-surfaces at time step 120,000; a 
x — z cross section slice is also shown. 


5.1 Shock-Cell Structure 

Experimental results for jets are often documented in 
terms of Schlieren pictures. It is straightforward to con- 
struct Schlieren plots (density-gradient modulus) from 
the numerical results. For the case of Mj = 1.42, 
Fig. 7 shows the experimental [5, Fig. 4(a)] and nu- 
merical Schlieren pictures. Good agreement in shock- 
cell structure is shown. For example, the shock cell 
width (spacing) is about 1.2877 in the streamwise direc- 
tion. Note that the experimental Schlieren plot is a time- 
averaged results, while the numerical one is a snapshot at 
a relatively early stage (time step 1 10,000). In the exper- 
imental Schlieren plot, it is observed that the first shock 
cell appears to be sharp and clear since the shear-layer in- 
stability wave is too weak at this location to significantly 
affect the shock cell. However, once the instability wave 
has gained a sufficient amplitude through its streamwise 
growth, it interacts strongly with the shock cells. This 
is evident from the deformations from the second shock 
cell and onwards and the additional blurring due to time- 
averaging in top panel of Fig. 7. 

5.2 Near-Field Radiating Screech Waves 

Figs. 8- 12, which represent a series of instantaneous 
snapshots of pressure iso-surfaces in the flow field, illus- 
trate the generation and propagation of screech waves. 
Since no forcing is applied in the numerical simulation, 
these waves are a clear indication of a self-sustained os- 
cillation. For Mj = 1.42, the jet screech is in the flap- 
ping B mode, which is a truly three dimensional one. The 
screech waves not only propagate in the upstream direc- 
tion but also swirl and flap around the jet core. 


After the computation has run for about 1 10,000 time 
steps, truly three dimensional asymmetric wave patterns 
begin to appear. Fig. 8 and Fig. 9 demonstrate the pres- 
sure iso-surfaces at time steps 120,000 and 130,000 re- 
spectively and from different angle. The iso-surfaces 
clearly take on nonaxisymmetric forms. However, at 
time steps of 200,000 and 320,000, the pressure iso- 
surfaces again turn out to look more axisymmetric, as 
demonstrated in Fig. 10 and Fig. 11 respectively. At 
time step 510,000, when the computation was stopped, a 
highly asymmetric pattern can be observed, see Fig. 12. 

5.3 Screech Frequency and Sound Pressure 
Level (SPL) 

The numerical time history is recorded for the location 
(2. 0,0. 6,0) at the nozzle lip in the flow field and later 
post-processed to obtain spectral information using Fast 
Fourier Transform (FFT) techniques. The recording be- 
gins after an initial time period has elapsed (180, 000 ac- 
tual time steps) ensuring that any start-up transients have 
left the computational domain. Figure 13 displays the 
SPL for the case of Mj = 1.42. The SPL plot shows 
the spikes at the fundamental ‘B’ mode frequency and 
its harmonics. The fundamental frequency spike corre- 
spond to screech frequency of 4500 Hz (SPL=122 dB), 
which agrees quite well with the experimental value of 
4350 Hz. 

5.4 Performance of the Parallel Computation 

The computation was mainly carried out on a Linux Pen- 
tium III PC cluster (NASA Glenn Research Center CW-7 
cluster) using 60 processes, running on 20 CPUs. Nor- 
mally, one would want a one-to-one relationship between 
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Figure 9: Pressure iso-surfaces at time step 130,000. 



Figure 10: Pressure iso-surfaces at time step 200,000. 



isffsurface levels: .7132, .71355, .7137 


Figure 1 1 : Pressure iso-surfaces at time step 320,000; a 
x — z cross section slice is also shown. 


the number of processes and processors. This less than 
optimal situation lead to about 15 % increase in the wall- 
clock time and was due to the computation having been 
started on a different Linux PC cluster using 60 CPUs. 
In the optimal situation, it takes about 8 seconds wall- 
clock time to march one step on the CW-7 cluster. The 
parallel computation was also tested on the SGI (Sili- 
con Graphics ) Origin 3800 workstation cluster with 64 
and 128 processors at NASA Ames Research Center. 
From 64 to 128 processors, the clock time reduces lin- 
early with the increasing number of processors, but the 
reduction will begin to level off if more processors are 
used. As a result of the explicit time-marching in the 
scheme, when running with 64 processors, the number 
of MFLOPS (megaflops) per processor remains between 
170 and 183, exceeding 20 % (160) of the theoretical 
peak MFLOPS. This performance is considered excel- 
lent by code-performance specialists at NASA Ames. 

6 Concluding Remarks 

In this paper, a 3-D CE/SE N-S solver using an un- 
structured hexahedral grid is briefly described and tested 
in a 3-D circular jet screech problem. The use of a hexa- 
hedral grid rather than a tetrahedral one enhances the nu- 
merical stability of the scheme and significantly reduces 
both memory size and CPU time, making the CE/SE 
scheme a viable tool for near-field CAA simulation. 

For the test case of M 3 = 1.42, the jet screech fre- 
quency of the dominant nonaxisymmetric ‘B’ mode and 
the shock cell structure agree well with the experimen- 
tal data [5, 6], Perhaps due to the relatively coarse grid 
used, the computed SPL is somewhat lower than the ex- 
perimental one. Further tests with refined grids are being 
carried out and will be reported in the future. 
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Figure 12: Pressure iso-surfaces at time step 510,000 
viewed from two different angles; x — z cross section 
slices are also shown. 


CE/SE simulation: MJ = 1 .49 — simplified nozzle geometry, 3.7M cells 
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